#include "navierstokes.hpp"

#include <functionSup.cpp>

#include <ns_cl.cpp>

namespace Feel
{
namespace FSI
{

template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
void
NavierStokes<Dim, Order, Entity>::updateResidual( const vector_ptrtype& X, vector_ptrtype& R )
{

    using namespace Feel::vf;

    if (M_LookComment) std::cout <<"[NavierStokes] : updateResidual start \n";

    boost::timer time;
    time.restart();

    this->changeRepository();

    //--------------------------------------------------------------------------------------------------//

    auto mesh = this->mesh();
    auto Xh = this->functionSpace();

    element_fluid_type U( Xh, "u" ); U = *X;
    element_fluid_type V( Xh, "v" );
    element_fluid_0_type u = U.template element<0>();
    element_fluid_0_type v = V.template element<0>();
    element_fluid_1_type p = U.template element<1>();
    element_fluid_1_type q = V.template element<1>();

    //--------------------------------------------------------------------------------------------------//

    double EE = 21*1e5;
    double sigmaS = 0.28;
    //double mu = EE/(2*(1+sigmaS));
    double muS = EE/(2*(1+sigmaS));
    double lambdaS = EE*sigmaS/((1+sigmaS)*(1-2*sigmaS));
#if 0
    auto Iddd = oneX()*trans(oneX()) + oneY()*trans(oneY());
    auto Fv = Iddd + gradv(*M_aux);
    auto Ev = 0.5*(gradv(*M_aux)+trans(gradv(*M_aux)) ) + 0.5*trans(gradv(*M_aux))*gradv(*M_aux);
    auto Sv = lambdaS*trace(Ev)*Iddd + 2*muS*Ev;

#else
    auto Iddd = oneX()*trans(oneX()) + oneY()*trans(oneY());
    auto epsv = 0.5*(gradv(*M_aux)+trans(gradv(*M_aux)));
    auto Sv = lambdaS*trace(epsv)*Iddd + 2*muS*epsv;
#endif

    //--------------------------------------------------------------------------------------------------//

    //Tenseur des deformations
    auto defv = 0.5*(gradv(u)+trans(gradv(u)) );
    auto def = 0.5*(grad(v)+trans(grad(v)) );
    //Matrice Identitee
    auto Id = oneX()*trans(oneX()) + oneY()*trans(oneY());
    // Tenseur des contraintes (trial)
    auto Sigmav = -idv(p)*Id + 2*M_mu*defv;
    // Tenseur des contraintes (test)
    auto Sigma = -id(q)*Id + 2*M_mu*def;
    //force volumique :
    //auto f = -rho*cst(9.81)*vec(cst(0.),cst(1.));
    auto f = vec(cst(0.),cst(0.));

    //--------------------------------------------------------------------------------------------------//

    form1( Xh, R, _init=true );

    bool M_isNLTerme=true;
    if (M_isNLTerme)
        if (M_isMoveDomain)
            {
                //le terme non linéaire
                form1( Xh, R ) +=
                    integrate( elements(mesh), _Q<3*(Order)>(),
                               //rho*trans(idv(u))*gradv(u)*id(v) );
                               //rho*trans(idv(u) -idv(*M_mesh_velocity))*(gradv(u))*id(v) );
                               rho*trans(id(v))*(gradv(u))*(idv(u) - idv(*M_mesh_velocity)) );
            }
        else
            {
                //le terme non linéaire
                form1( Xh, R ) +=
                    integrate( elements(mesh), _Q<3*(Order)>(),
                               rho*trans(id(v))*(gradv(u))*(idv(u) ) );
            }
#if 1
    // sigma : grad(v) sur Omega
    form1( Xh, R ) +=
        integrate( elements(mesh), _Q<2*(Order)>(),
                   + trace(Sigmav*trans(grad(v))) );
#if 0
    form1( Xh, R ) +=
        integrate( elements(mesh),// _Q<2*(Order)>(),
                   - rho*divv(*M_mesh_velocity)*trans(idv(u))*id(v) );
#endif
#else
    form1( Xh, R ) +=
        integrate( elements(mesh),
                   + 2*M_mu*trace(defv*trans(grad(v))) );
    form1( Xh, R ) +=
        integrate( elements(mesh),
                   + div(v)*idv(p) );


#endif

    // terme assurant l'incompressibilité
    form1( Xh, R ) +=
        integrate( elements(mesh), _Q<2*(Order-1)>(),
            + divv(u)*id(q) );

    auto clDef = CL_DEF_NS(M_time);
    ForEachCL( clDef,cl::neumann_scal,
               form1( Xh, R ) +=
               integrate( markedfaces(mesh,mesh->markerName(PhysicalName)),
                          - trans (Expression*N())*id(v) )     );

   //--------------------------------------------------------------------------------------------------//
    if (M_weakCL)
        {

            ForEachCL( clDef, cl::dirichlet_vec,
                       form1( Xh, R ) +=
                       integrate( markedfaces(mesh,mesh->markerName(PhysicalName)),// _Q<2*Order>(),
                                  - trans(Sigmav*N())*id(v)
                                  + penalbc*trans(idv(u)/*-Expression*/)*id(v)/hFace() ) );

            ForEachCL( clDef, cl::paroi_mobile,
                       form1( Xh, R ) +=
                       integrate( markedfaces(mesh,mesh->markerName(PhysicalName)),// _Q<2*Order>(),
                                  - trans(Sigmav*N())*id(v)
                                  + penalbc*trans(idv(u)/*-idv(this->meshVelocity())*/)*id(v)/hFace() ) );

#if 0
            // terme de bord classique
            form1( Xh, R ) +=
                integrate( markedfaces(mesh,mesh->markerName("Entree")), _Q<2*Order>(),
                           - trans(Sigmav*N())*id(v) );

            if(M_isFixeBoundary)
                form1( Xh, R ) +=
                    integrate( markedfaces(mesh,mesh->markerName("Paroi")), _Q<2*Order>(),
                               - trans(Sigmav*N())*id(v) );

            // Vitesse d'entree du fluide
            auto Vit_Entree = 10*Vitesse_Entree*( Py()*(1.0-Py()) )*oneX();
            //terme de bord weak boundary
            form1( Xh, R ) +=
                integrate( markedfaces(mesh,mesh->markerName("Entree")), _Q<2*Order>(),
                           - trans(Sigma*N())*(idv(u)-Vit_Entree)
                           + penalbc*trans(idv(u)-Vit_Entree)*id(v)/hFace() );
            if(M_isFixeBoundary)
                form1( Xh, R ) +=
                    integrate( markedfaces(mesh,mesh->markerName("Paroi")), _Q<2*Order>(),
                               - trans(Sigma*N())*idv(u)
                               + penalbc*trans(idv(u))*id(v)/hFace() );
            form1( Xh, R ) +=
                integrate( markedfaces(mesh,mesh->markerName("ParoiH")), _Q<2*Order>(),
                           //- trans(Sigma*N())*idv(u)
                           + penalbc*trans(idv(u)-idv(*M_mesh_velocity))*id(v)/hFace() );

            form1( Xh, R ) +=
                integrate( markedfaces(mesh,mesh->markerName("ParoiB")), _Q<2*Order>(),
                           //- trans(Sigma*N())*idv(u)
                           + penalbc*trans(idv(u)-idv(*M_mesh_velocity))*id(v)/hFace() );
#endif

        }
   //--------------------------------------------------------------------------------------------------//

    //Force Volumique
    form1( Xh, R ) +=
        integrate( elements(mesh), _Q<Order>(),
            - trans(f)*id(v) );


    //bool M_isMoveDomain=false;
    if (M_isMoveDomain)
        {
#if 0
            form1( Xh, R ) +=
                integrate( markedfaces(mesh,mesh->markerName("ParoiH")), _Q<3*Order>(),
                           - trans(Fv*Sv*N())*id(v) );
            form1( Xh, R ) +=
                integrate( markedfaces(mesh,mesh->markerName("ParoiB")), _Q<3*Order>(),
                           - trans(Fv*Sv*N())*id(v) );
#else
            /*//new comment ??
            form1( Xh, R ) +=
                integrate( markedfaces(mesh,mesh->markerName("ParoiH")), _Q<3*Order>(),
                           - trans(Sv*N())*id(v) );
            form1( Xh, R ) +=
                integrate( markedfaces(mesh,mesh->markerName("ParoiB")), _Q<3*Order>(),
                - trans(Sv*N())*id(v) );*/
#endif
        }
    else if (M_weakCL)
        {
#if 0
            // terme de bord classique
            form1( Xh, R ) +=
                integrate( markedfaces(mesh,mesh->markerName("ParoiH")), _Q<2*Order>(),
                           - trans(Sigmav*N())*id(v)
                           + penalbc*trans(idv(u))*id(v)/hFace() );
            // terme de bord classique
            form1( Xh, R ) +=
                integrate( markedfaces(mesh,mesh->markerName("ParoiB")), _Q<2*Order>(),
                           - trans(Sigmav*N())*id(v)
                           + penalbc*trans(idv(u))*id(v)/hFace() );
#endif
        }
   //--------------------------------------------------------------------------------------------------//

    /*    if (M_weakCL)
        {

            //terme de bord weak boundary
            form1( Xh, R ) +=
                integrate( markedfaces(mesh,mesh->markerName("Entree")), _Q<2*Order>(),
                           + trans(Sigma*N())*Vit_Entree
                           - penalbc*trans(Vit_Entree)*id(v)/hFace() );
            form1( Xh, R ) +=
                integrate( markedfaces(mesh,mesh->markerName("ParoiH")), _Q<2*Order>(),
                           //+ trans(Sigma*N())*idv(*M_mesh_velocity)
                           - penalbc*trans(idv(*M_mesh_velocity))*id(v)/hFace() );
            form1( Xh, R ) +=
                integrate( markedfaces(mesh,mesh->markerName("ParoiB")), _Q<2*Order>(),
                           //+ trans(Sigma*N())*idv(*M_mesh_velocity)
                           - penalbc*trans(idv(*M_mesh_velocity))*id(v)/hFace() );
                           }*/
    //------------------------------------------------------------------------------------//

    //bool M_isStat=true;
    if (!M_isStat)
        {

            form1( Xh, R ) +=
                integrate( elements(mesh), _Q< 2*Order >(),
                           + rho*trans(idv(u))*id(v)*M_bdf_fluid->polyDerivCoefficient(0) );

            element_fluid_type Buzz( Xh, "v" );
            element_fluid_0_type buzz = Buzz.template element<0>();
            buzz.container() = M_bdf_fluid->polyDeriv().template element<0>().container();

            //auto buzz = M_bdf_fluid->polyDeriv().template element<0>();

            //auto buzz = M_bdf_fluid->unknown(0).template element<0>();
            //auto buzz = M_bdf_fluid->poly().template element<0>();
#if 0
            form1( Xh, R ) +=
                integrate( elements(mesh), _Q< 2*Order >(),
                           - rho*trans(idv(buzz))*id(v)*M_bdf_fluid->polyDerivCoefficient(0) );
#else
            form1( Xh, R ) +=
                integrate( elements(mesh), _Q< 2*Order >(),
                           - rho*trans(idv(buzz))*id(v) );
#endif

        }

    //------------------------------------------------------------------------------------//

#if 0
 int stabilisation = this->application()->vm()["stab"].template as<int>();
if (stabilisation==1) {
    //terme de stabilisation
    form1( Xh, R ) +=
        integrate( internalfaces(mesh), _Q< 3*Order >(),
           // + 0.025*vf::pow(hFace(),2.0)*( abs( trans(idv(u))*N() ))*(trans(jumpv(gradv(u)))*jump(grad(v)) )  );
            + (0.05+vf::pow(hFace(),2.0))*( (trans(idv(u))*N())*( trans(idv(u))*N() ))*(trans(jumpv(gradv(u)))*jump(grad(v)) )  );
}
#endif

    //------------------------------------------------------------------------------------//

    R->close();


    if (true)//!M_weakCL)
        {
            auto dirElem = this->getDirichletElement();
            ForEachCL( clDef,cl::dirichlet_vec,
                       modifVec222(*dirElem,R, vf::vec(vf::cst(0.),vf::cst(0.)) , PhysicalName ) );
            ForEachCL( clDef,cl::paroi_mobile,
                       modifVec222(*dirElem,R, vf::vec(vf::cst(0.),vf::cst(0.)) , PhysicalName ) );
        }


    double ti = time.elapsed();
    if (M_LookComment) std::cout <<"[NavierStokes] : updateResidual finish in "<<ti<<"s\n";



}

} // end namespace NavierStokes
} // end namespaceFeel


// instantiation
template class Feel::FSI::NavierStokes<2,2,Feel::Simplex>;
//template class Feel::NavierStokes<2,4,Feel::Simplex>;
